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Abstract 

Conformal arrays are popular antennas for aircraft and missile 
platforms due to their inherent low weight and drag properties. How- 
ever, to date there has been a dearth of rigorous analytical or nu- 
merical solutions to aid the designer. In fact, it has been common 
practice to use limited measurements and planar approximations in 
designing such non-planar antennas. In this report, we extend the fi- 
nite element-boundary integral method to scattering and radiation by 
cavity-backed structures in an infinite, metallic cylinder. In particular, 
we discuss the formulation specifics such as weight functions, dyadic 
Green’s function, implementation details and particular difficulties in- 
herent to cylindrical structures. Special care is taken to ensure that 
the resulting computer program has low memory demand and mini- 
mal computational requirements. In this report, both scattering and 
radiation parameters are computed and validated as much as possible. 
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1 Introduction 


Conformal antenna arrays are attractive for aircraft, spacecraft, and land ve- 
hicle applications since these systems possess low weight, flexibility, and cost 
advantages over conventional protruding antennas. The majoritv of previous 
studies in non-planar conformal antennas has been conducted experimentally 
due to a lack of rigorous analysis techniques. Some approximate analyses 
have been considered but these are restricted in accuracy and element shape. 

Recently, the finite element-boundary integral (FEM-BI) method was 
successfully employed for the analysis of large cavity-backed planar arrays 
[1]. The resulting system is sparse due to the local nature of the finite ele- 
ment method whereas the boundary integral sub-matrix is fully populated. 
However, by resorting to an iterative solver such as the Biconjugate Gradi- 
ent (BiCG) method, the boundary integral system may be cast in circulant 
form allowing the use of the Fast Fourier Transform (FFT) in performing 
the matrix- vector products. This BiCG-FFT solution scheme ensures C?(N) 
memory demand for the entire FEM-BI system and minimizes the computa- 
tional requirements. 

In this report we extend the FEM-BI formulation to aperture antennas 
conformal to a cylindrical metallic surface. Both the radiation and scatter- 
ing problems will be developed in the context of the FEM-BI method. In 
contrast to the planar aperture array, the implementation of the cylindrically 
conformal array requires shell-shaped elements rather than bricks, and the 
required external Green's function must satisfy the boundary conditions on 
the cylinder. In its exact form, this Green’s function is an infinite series which 
imposes unacceptable computational burden on the method. However, for 
large radius cylinders, suitable asymptotic formula are available and herein 
used for an efficient evaluation of the Green’s function. In addition, the 
resulting BI system is again cast in circulant form to ensure an O(N) mem- 
ory demand and take advantage of the FFT’s speed when carrying out the 
matrix-vector product. 

At the end of this report, we also present computations of the relevant 
antenna and scattering parameters validated with available known results. 
One of the most interesting applications of the numerical laboratory devel- 
oped in this project will be the ability to study curvature effects of real world 
antenna systems in an exact manner. Such studies will be pursued in the 
near future and the will be addressed in a future report. 
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2 FEM-BI for Circular Cylinders 

Consider a cavity recessed in an infinite.metallic cylinder, shown in figure 1. 
The cavitv walls are assumed to coincide with constant p-, Q- and z-surfaces 
and the cavity is assumed to be filled with inhomogeneous material, nter.or 
sources and lumped loads may be present as well as surface metallization 

patterns and resistive cards. _ . . . , , 

The FEM-BI technique [1] permits determination of the electric fi 

within the cavity which are induced by either interior or exterior sources. The 
svstem of equations associated with the FEM formulation is sparse whereas 
the boundary integral sub-matrix is usually fully populated However, wit 
a judicious choice of boundary elements, the formulation will maintain O(N) 
memory demand when coupled with a Biconjugate Gr ^ ,e “ t *F ast 7°^ 
Transform (BiCG-FFT) solver. Upon determination of the fields, the radi- 
ated and scattered patterns are readily calculated from the aperture fields 
while the input impedance. S-parameters and other pertinent antenna quan- 
tities mav be computed from the appropriate interior fields. 

The FEM-BI system is developed directly from the inhomogeneous vec or 
wave equation and a complete presentation of the derivation is given y 
Volakis et al. [2]. The resulting system is 

r $ V x Wj(p, <fr, ;) • V x Wj{p, <P, z) 

JV, | pr{Pi4 > -< z ') 

4 >,z)-W i (p,<i>,z)\pdpi l t>iz 

+(* ,a)X(j)6.(<) J $ J 

Zi 2 (a,<p,z) x /3(a, <j>',z) ■ Wj{a,<f>',z )] d<j) dz d<t>dz = f\ nt + /, (1) 


where tU,(p,d>, r) are edge-based expansion functions where support is over 
the volume V, and £ 2 (a,<M) is the pertinent dyadic Green’s function. The 
free-space wavenumber is denoted by K = % the cylinder has radius a 
whereas e r and p T are the relative constitutive parameters of the material 
filled cavitv. Also Vj represents the i th integration volume which extends 
over the support of VF,(p,<M) and in a similar fashion 5, and S, indicate 
surface integration over all aperture elements which extend over the support 
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of the i th (test) and j th source basis functions, respectively. The function 
6,,{i)fi a (j) is a product of two Kronecker (unctions and it simply indicates 
that the boundary integral only contributes when both the test and source 
unknowns are on the aperture. We remark that the dyadic Green s function 
is convolutional (0 = 0 - 0 '.: = =-=') and will be discussed later along 
with the functionals f\ nt and ff zt which depend on the interior and exterior 
excitations, respectively. Rather, we shall first proceed with a presentation 
of the appropriate weight functions and the evaluation of the FEM matrix 
entries which are represented by A,-j in the FEM-BI matrix s\stem 


A] 



+ 


IS] (01 
[ 0 ] [ 0 ] 


{cr} _ ur ) 


( 2 ) 


to be solved via the BiCG method. 


2.1 Vector Weight Functions 

Traditional node-based finite elements associate the degrees of freedom (£,) 
with the nodal fields. These elements impose both tangential and normal 
field continuity along inter-element boundaries even at material boundaries 
[3], These elements also do not correctly model the null space of the curl 
operator and as a result spurious solutions are generated [3] which are typ- 
ically suppressed by imposing a penalty function [4]. Edge-based elements, 
where the unknowns are associated with the free edges of the mesh, do not 
enforce normal field continuity (although tangential continuity is still main- 
tained) and are therefore more suitable for electromagnetics applications. In 
addition, edge-based elements avoid explicit specification of the fields at ge- 
ometry corners where the field may be singular [5]. Two popular volume 
elements are the brick [5] and the tetrahedron [6, 7] which are readily mated 
with rectangular and triangular surface elements, respectively. The brick is 
well-suited for geometries delimited by constant x- , y- and z-planes in the 
Cartesian system while the tetrahedra are more versatile (and consequently 
more complex). For cavities residing in a circular cylinder it is advantageous 
to employ the cylindrical shell element. 

Cylindrical shell elements, which are defined by constant p-, <f>- and z- 
surfaces, are superbly suited for cavities recessed in a circular cylinder since 
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t hey possess both geometrical fidelity as well as simplicity. These shell ele- 
ments are analogous to the bricks used by Jin and Volakis [oj and belong in 
the general class of curvilinear brick elements. Figure 2 illustrates a typical 
shell element which has eight nodes and twelve edges: four edges aligned 
along each of the three orthogonal directions of the cylindrical coordinate 
system. The edge-based weight functions for these elements must satisfy the 
following properties: 

1. subdomain (finite element) 

2. || o , :) ||= 0 if (p,0,z) € any edge ^ j and || j ,h edge 

3. V - 1 Vj{p.o.z) = 0 1 

Such expansion/ weight functions can be represented as 

vT'i 2 (p.6,z) = \X'„(p. 0 .Z;-. 0 r ,Z t .+), W<z{p. 0 ,z) = w„(p. 0 ,ZX-, 0 l,Zt,-) 
TVs e(p,0,-) = lTp(/3, <?, z;-<<pr,Zb , -), -) = W p (p,<p,z;-,<f>i,zt > ,+) 


W n( p, s ) — TV 0 ( p- - : pb y ■ ■ ~t ' T ) i TV 23 ( Pi 0i ~ ) TV ^,( p, 0, ~ , p a , • , ) 

0. -) = VT«>(p, 0 . c; pb, •, Zb, — ), W 6 7{p,4>,~) = VF^p, <j>, z\ p a , •, Zb, + ) 


^Vis(p, C) — TV - ( p, 0, , P(», (^ r , ’ , T VV26(/b ) VV 2 (p, 0, Pa, <pr, •, ) 

TT 48 (p, o, r) = VF.(p, c£), r; pt, <pi, •, -), VV^p, 0 , 2 ) - W z {p, 0, z\ p a , 0/, •, + ) 


(3) 


where TT^. denote the edge which is defined by local nodes (l,k) as shown 
in figure 2. That is three fundamental vector functions are required for the 
complete representation of the edge-based cylindrical shell element. They 
are given by 


W p {p,4>,z-,p,0,z,s) 


sfb (0 ~ 0)U ~ =) - 

ah p 


1 \V } ( p , 0 , c) will only satisfy this requirement within the volume of the element. These 
weighting functions introduce artificial charges on the faces of the element and are not 
divergenceless at element interfaces. This is required since these elements do not guarantee 
normal field continuity across the element faces. 
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(- 1 ) 


11 «>(/>. o. z: p, o. i, i) 
U'-(p, o, z; p , o, z, 5) 


-r)o 

r-(/ J - ^)(o - 


u here t P6 p a , q — c> r o/, and /i — — z^ and the element parameters 

(pa- Pb,<PhO r . Zb. z t ) are shown in figure 2. We remark that the - term which 
appears in the definition of the p-directed weight is essential in satisfying the 
divergenceless property of the edge-based expansions. Of course, for very 
large radius cylinders and small elements, the curvature of these cylindrical 
shell elements decreases resulting in weight vectors which are functionally 
similar to the bricks used by .Jin and Volakis [5]. 


2.2 FEM Matrix 


The inherent locality of a partial differential equation formulation such as 
the FEM suggests that [A] is sparse. This FEM matrix results from the 
discretization of the first two integrals of (1) and owing to the finite support 
of the basis functions its entries are identically zero unless both the test and 
source edges are within an element. These matrix entries can be represented 


as 


4 — _ r( , )b 7,2 /( 2).j 

‘ M.7 — *st ~ 

Pr 


(5) 


where we have assumed constant material properties within an element (e r 
and p T ) and the subscripts (i,j) refer to the row and column of the [A] matrix, 
respectively. The two auxiliary functions are defined by 


- f V X W,(p , <p, fa, (pj, Zj,s } ) ■ 

J V , 

V x W t (p,<j>,z’, pi,<t>i,z t ,Si)pdpd<f>dz 
= J y W a {p,<i>, p } ,<}>j, Zj,Sj) ■ W t (p, <p, z-,p t , <t>i,z t ,s,)pdpd(pdz ( 6 ) 


where (s,t) € {p,<j>, z} indicate the direction of the source and test edges, 
respectively. Since the fundamental weight vectors (4) are aligned along 
orthogonal directions, (6) is symmetric with respect to the source and test 
edges (e.g. l[ t * 3 = / and thus, only six combinations of (s,t) for 1^ 
need be determined and only three such combinations are required for I^\ 
They are 
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i 


/(>) = 
pp {ah ) 2 


p\h In 



Pb 

. Pa J ^$1 


<Pr 


(0 - <£,)(<? - 0*)d0 + 





(*-£.)(*-■ 


/<!> = - 1 


^s$tPb 


■ f {4> - O s ){<? - <Pi)dQ 
a i J<pi 

(2 {A - o2) - » (ft + ft) + ftft ln (*)) /* ( * " 

il\ ] = j P \p - ps){p - pMp 

e - ^[«(5W-^)-“A +w+ '* Aln (* 

5 (pi - pi) J* w - 

r(2) = r (4> - isM - 4>t)d4> [ ‘ {Z -*•)(■■ 

pp {ah ) 2 \Pa) hi Jzb 

1 % - (^ -/>:)+ 5 ( ft+ft>('’--^) + 5^‘ 

/*'(* - *,)(« - 5|)<k 

/«> = ^[j((>!-(-l) + 5(ft + ft)('’--' , ‘) + 5^' 

f {<b — 4> s ){<p — 4>t)d4> 

J<h 

where each of the unevaluated integrals are of the form 

r ({ - L) (f - «<) <*e = | ( i2 - (,J ) («• + *■) + 3 ~ 

J Ij 


z t )dz 

Pi ( pi - pi}) + 

z a )(z-z t )dz 

+ 

z — Zf)dz 

{pi-pi)] x 
W-P 2 a )] x 

(7) 

I 3 ) + tit (U - L) 

( 8 ) 
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I lie [A] matrix is assembled by evaluating (5) for each element combination 
which contain the i th and j' h unknowns. Since the integrals (6) are svm- 
iih tiic. onh the upper triangle or lower triangle of the FEM matrix need lie 
stored. 


2.3 Boundary Integral Matrix 

The entries of the boundary integral sub-matrix, [(?], couple the interior 
fields of the FEM system and any external impressed fields. In addition, the 
boundary integral provides an exact boundary condition for mesh closure. 
The entries of [Q] are 


Gtj = (k Q a ) 2 J s J s lT,(a.C>. z:p,. o,, i,, i.) • 

p(a,o,:) x G’ 2 (a,0, =) x p(a,<j>',z) 

•11 4 (n. 0 , z ; p-,, pj, ij, Sj)d<t> dz'dcpdz (9) 

and we note that the global nature of the Green’s function implies a fully 
populated sub-matrix (e.g. all boundary unknowns couple to all the other 
boundary unknowns). The surface weight functions in (9) are the volume 
weight functions (3) collapsed down to a cylindrical-rectangular patch on the 
surface of the cylinder, p = a. The dyadic Green’s function ^ 2 ( 0 ,^, z) is of 
the second kind and it satisfies the Neumann boundary condition 


p x V x G 2 {a,4>,z) = 0. (10) 

on the cylinder s surface as well as the Sommerfeld radiation condition. It 
can be expressed component-wise as 


G 2 {a, 0 ,z) = ii'G*(a,hg)+fe + z$]G+*(a,4>,2) + zzG M (a,faz) 

( 11 ) 


where the unprimed unit vectors are functions of the test point (a,^,z) and 
the primed unit vectors are functions of the source point (a y <f>\z). The 
components may be exactly expressed as an angular eigenvalue series [8] 


G zz (aJ,z) 


1 

(2a 0 2 



1 ^ 2) ( 7 ) 

-< hT\ 7 ) 


e 


j(ni ~ k ‘ z) dk 2 
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G°*{a,o.:) 


1 V~' r ( nkz \ 

C~la.o.:) - -JZyiJ^L *?>(,) 

* rcc 1 

S.L; 


C2-) 2 


HT\ 7) ( nk, V W h 


//i 2 , b) v Av/lA ' 


//n 2 V; 


e j(n *~ k ‘ i) dk z 


(1 


However, for large radius cylinders, (12) is computationally prohibitive. 
Instead, it is advantageous to employ one of the available asymptotic approx- 
imations ofS,(a.*,5) [9, 10. 11, 12]. These approximations improve as the 
cylinder radius and the geodesic distance between the source and observation 
point increase. Since each approximation invokes Watson s transformation, 
thus converting the angular eigenvalue series of (12) into a creeping wave se- 
ries, the geodesic or on-surface ray parameters are necessary in des ««bmg e 
interaction between surface fields. The formula due to Pathak and Wang [9] 
have proven most accurate for our applications since the regions where Bird s 
formula [12] and Boersma and Lee’s expressions [10] are superior (paraxia 
and far-zone) correspond to situations where coupling is small. In add, ion, 
we shall find it necessary to employ a regularization process to correct for 
errors in computing the near region (s - 0) and to also allow evaluation 
of the singular integrand as (<p,r) -(*,*) which mitigates any advantage 
t hat the more complex formula have with respect to Pathak s expressions in 
the near-zone. Using a uniform expansion of the Hankel functions in ( ) 

and a first-order evaluation of the axial wavenumber integral, Pathak found 

that v 

[)(2 — 3cos 2 #)) v(0) i 

1 


G’”(a,<?, z) ~ 


(cos 2 0 + q( 1 

G* z (a,i>,~) ~ 

i^-qe~ ]ko3 sinflcosfl j (1 — 

G**{a,4>,z) ~ 

-itl qe -i k ° 3 < 

2tt q 



+q ^sec 2 # {u(0) — u(/3))] | 


( 13 ) 


where 0 = ks [$&] 5 and tlie g eodesic parameters are given by 
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Ceodrsi 


>ic 


path length \ s = yf(a&) + -’j 


2. Geodesic trajectory (o = tan~ l |is±lj j 

3. Approximation order (r/ = 

where $ = o or $ = 2 jt - q depending on which of the two direct paths 
illustrated in figure 3 (the short or long) is used. In (13). u{3) and v™(3) 
represent the soft and hard surface Fock functions, respectively. These func- 
tions are characteristic of the creeping waves on a circular cylinder and are 
discussed in detail by Logan [14] along with their evaluation. 

As stated previously, a regularization process is necessarv for the accurate 
evaluation of (9) especially when the source point and test point coalesce. 
Bird has shown [12] that Pathak's approximation (13) reduces to the metallic 
screen Green s function. 2G 0 (a, <2>. ~), modulated by v{0) within the available 
order (D{q )). This suggests a regularization procedure where the metallic 
screen Green's function is subtracted from the cylindrical Green's function, 
thus rendering it non-singular, followed by an addition of the planar contri- 
bution. The non-singular Green's function is given by 

G~ z (a, q, ~) ~ 


G <t>z (a , 0, c ) ~ 


<p, c) ~ 


(14) 

and this asymptotic form of the Green’s function is used for all short path 
contributions whereas (13) is used for the long path contribution in which 
case curvature effects are expected to be prominent. When the regularized 
integrand is used, the contribution of the planar Green’s function 


~ J l^ qt Jk °\ ( C0S *° + - 9)(2 - 3cos 2 0)) [t>(0) - 1] 

^e-^sinflcostfjo ~3r/(l -g))[i>(/?)- 1] j 

~T^ (ie ^ j {sin 2 0 + q{\ - q)( 2 - 3 sin 2 0)^ [v((3) - 1 ] 
+<7 [a ec 2 0(u(/3) - v(/3))j | 


G 


*J 


- 2{k 0 a) 2 j j W t (a, <p, c; />,, i, ) • 
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p(o. o. z ) x 0 o (a. o. r ) x p(a, 0 , z ) 


where 


is added to (14) 


G 0 (a. o, c) = 


7 

k 2 

A o j 


Gij — G,j + G'f 


)do dz dodz 

(15) 



4 -R 

(16) 


(17) 

ix + yy + zz. 

Upon use of a 


(15) 


G l} 


( k 0 a 


'° a > f f i7- - = 

2 *~ Js, Js } n ^ a ' 0, P " ^ ' 1 

•M ,(a. o ,z :pj,<?j, Zj, Sj)dO dz do dz 

~2w IJs V - U 0,0, x Wi(a,0,c;/>,,0, ,£,,£,)] 


V • p(a , 0 , z') x W a ( a , <p'. 


1 

Pj • ■*> )J — — dtp dz d<p dz 

(18) 


and this form of the boundary integral may be readily evaluated even as R 
vanishes by employing the regularization procedure used by Jin and Volakis 

in C n< ? e that . 1 ’^ 1 as a oo and hence the regularized integrand 
(14) vanishes leaving only the planar contribution (18). With the specifica- 
t.on of both the FEM matrix. [A], and the boundary integral sub-matrix, 
[y\, the unknown fields can be determined for a given excitation. In the next 
section, explicit evaluations of the excitation functions appropriate for plane 
wave scattering and probe-fed antenna applications are given. 


3 Excitation 

Two sources of electromagnetic fields are considered: external sources (plane 
waves) for scattering analysis and internal sources (probe-feed) for antenna 
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parameter calculations. The use of the exact boundary condition in (1) 
allows the coupling of an exterior excitation field into the cavity while the 
FE.\1 formulation itself [2] readily permits modeling of "'I"'" 1 ' J" 

this section we describe the form of the source lunctionals /, and/, al g 

with their numerical implementation. 

The forcing functional. /“'* due to exterior sources is given by 


/?*‘ =jZ 0 k 0 aj ' s l .z)- p{a,Q <=) x H cyl {a,<P .z )do dz 


(19) 


where lT’(,i.o,r) is the testing vector for the i<‘ row of the matrix and it* 
represents is the magnetic held on the cylinder’s surface m the absence of 
the cavity. For plane wave incidence 

c* c~ 

fr = Y 0 (l‘ X 

, - * , . n 1 J ko[p sin 8 X cos {<!>-<Pi) + z cosOt] 

_ y A' sin 7 cos 0, - 0 cos 7 - - sin 7 sin 0 , j e 

‘ 1 ( 20 ) 

where (0„o,) indicate the direction of incidence, 7 is the polarization angle 
and e‘ = 0’cos7 + £sin7. The total surface field is given by 

( 21 ) 


H cyl (a. <p. z) = H‘{a,(p,z) + H* yl {a,<j>,z) 


with 


/ \ / sin 7 jfc 0 cosd,s \ v 


gj n ( f + 4 > ”^ > * ) 


jk Q cos 0 t z c* 


Hf(a,<p,=) = -2V; 

J 


nk 0 a sin 0, 
n sin 7 cos 0, 


n =-oo lf/n 2) (fcoasin0,) 
cos 7 


L 


li/ri 2) (A: 0 asin^) 
jn(* ) 


A’ 0 a sin Hn 2 ) (k 0 a sin 0,) J 


(22) 


These expressions may be computed by summing only a few terms of the 
series if k 0 a sin0, is small. However, as this parameter becomes large (e.g. 
for large a and 0, - 90°), (22) may be replaced with equivalent asymp- 
totic representations similar to those considered earlier. Utilizing Watson s 
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transformation and Fock theory [13] in connection with (22)* we find that 

2 

//?'’ ~ -Y 0 sin -y sin 0 l t jk ° cos6 ' : £ e -J^» [ 5 (0) ( W <I) ;> | 


p=i 


// 


cy/ 


j’2Vo cos 7 - — V* c J fc ° asine ' < * , p [/ ,0) ( /??$>„)]’ 
A*o«sin0, ^ L P/ J 


-Y 0 sin 7 cos 0,e jk ° cos6 ' z ^(-1 

p=l 


5 (O) (77?0 p ) 


/7? 


-77 — ^r</ (1, (”^ 


A\,fl sin 0 , 


(23) 


in which <fq = ^ - (o - o,), 4 > 2 = (p - 0 .) 


m = 


' , and the 


symbol indicates complex conjugation. The appropriate Fock functions are 

[M ] 2 


9 W U) = 


f 


(0 


J— [ — 

v/T Jr w\ 

(0 = 4= J 

\J 7t 


(0 




«’i(0 




(24) 


where u'i(t) and its derivative tCj(t) denote the Airy functions and the inte- 
gration contour is given by Logan [14]. 

The asymptotic formulas (23) are quite accurate except when <t> ~ <p t . In 
this region, Goriainov’s [15] expressions 

II c . yl ~ — sin asin 0,e jfcoCO * fl ' : /e _ - ,/ ‘°" sl,10,<I>1 [< 7 * 0 *(m<I>i )j 

jk 0 a sin 0, cos (<p—4 , ) CQg (<p _ | 

™ jk o cos0,*( e -jk o a,ine.* t )] * 

[ L J 

+ e jk oa s> n*. cos (*-*,) [F(_ mcos( ^_ | 


+e J 


H C 0 l ~ j-Yo COS Q 


“Logan [14] uses thee ^ * time dependency in the definition of these functions requiring 
the complex conjugation in (23) 
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( 25 ) 


with [14] 

G (/) ( 0 = <7 (,) (0^ 

F li) (0 = / (,, (f)e j ^ ( 26 ) 


have been found to be more accurate and can be used instead of (23). 

These surface field expressions may now be used to efficiently calculate 
the entries of the column vector {/,"'} via a numerical evaluation of (19). In 
particular, the modal series (22) is used when k 0 a sin 0, < 10 and either (23) 
or (25) for k o as\n0, > 10 as appropriate. 

Conformal antenna patches are typically fed by a microstrip line printed 
along with the radiator on the surface of the substrate The microstrip lines 
are in turn fed by a coaxial probe which originates behind the cavity as shown 
in figure 4. The patch may also be fed directly by the coax feed or through 
some form of aperture coupling. Nevertheless, the principal excitation for 
the system is given by 


fl 


mt 



A/ ,nt (p, 4>, z) 


+ jk 0 Z 0 J' nl {p 


>,<M) j 


Wi{p, <f>,z)pdpd<j> dz 


(27) 


where J int and A/ ,nt are the impressed electric or magnetic current densities 
representing the sources. For a radially (p) directed probe feed, the impressed 
monopole current located at (<p a ,z 3 ) is given by 
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(28) 


id 




which results in an excitation function (27) 
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( 29 ) 


With both types of excitation and the FEM-BI matrix now specified, the 
BiC’G method may be used to determine the unknown electric fields within 
the cavity. 


4 System Solution 

The FEM-BI system (2) may be solved using one of the popular direct meth- 
ods such as Gauss- Jordan elimination, Gaussian elimination and LU decom- 
position. Alternatively an iterative methods such as Hotelling s method, 
conjugate gradient method or the biconjugate gradient method may be used. 
We have chosen to use the symmetric form of the BiCG solver because it 
requires much less memory than a direct method and moie importantly it 
can be implemented in a manner which is computationally efficient (utilizing 
only one matrix-vector product per iteration). Although use of an iterative 
method such as the BiCG method can require more wallclock time than a 
direct solver when multiple right-hand sides are considered (such as is the 
case with backscatter calculations), memory demand was deemed the most 
critical and expensive resource. The BiCG algorithm is given in Appendix 
A and the efficient FFT-based calculation of the boundary integral matrix- 
vector product is discussed in Appendix B. In this section, we will present 
details of this iterative solver specific to this application. 

The BiCG algorithm requires one matrix-vector product per iteration as 
shown in Appendix A. This operation represents the bulk of the computa- 
tional demand of the method and requires 0{N 2 ) complex operations in the 
case of non-symmetric matrices. The matrix-vector product is carried out by 
executing the sum 

y[n] = [A] {x} = 51 A[n,n']x[?j'] n = 1,2,3, ..., N (30) 

n' = l 

and if the matrix is sparse, a storage scheme such as the Compressed Sparse 
Row (CSR) format may be used to reduce both the memory demand and 
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computational load. Lsing the C'SR scheme. (.JO) can be rewritten as 


7-[nJ 

i/[n] = [/l] {-r} = -4[e[n, 7 ? = 1,2.3 V (31) 

n' = l 

where r[n] is the number of non-zero entries per row of the matrix and e[n.r/] 
indicates which entry of the long data vector, .4. is associated with the matrix 
entry A [n, n ]. The FEM matrix, [A] in (2) is such a sparse matrix. Although 
additional memory savings are possible since the FEM matrix is symmetric 
(thus only upper triangle of the matrix need be stored), experience has shown 
that use of a symmetric matrix-vector product leads to a severe performance 
degradation on vector computers (such as a CRAY) due to the resulting short 
vector lengths. Therefore, the entire sparse FEM matrix is stored and used 
in the product so that the code s performance is maximized when executed 
on vector architectures. 

The boundary integral matrix-vector product involves the fully populated 
matrix, [£]. If uniform surface elements are used in the discretization, this 
matrix-vector product may be expressed as a truncated. discrete, linear con- 
volution and thus amenable to efficient calculation using the Fast Fourier 
Transform (FFT). Although uniform zoning imposes restrictions on the ge- 
ometries which can be analyzed by this FEM-BI technique, the resulting 
memory and computational efficiencies have proven to be well worth the 
saciifice. The boundary integral product is implemented as described in 
Appendix B with the following problem specific exception. The cross-term 
arrays do not possess the property: G*. [m - m', n - n'] = C^ z [m' - m, n - n] 
and hence the periodic replication rule (B-3S) cannot be used here. In- 
stead. a different replication rule was developed in performing the discrete 
convolution. The resulting code has proven robust even for large arrays. In 
addition, the cross-term replication rule must be changed for wraparound 
arrays since the cavity discretization is fundamentally different than is the 
case for cylindrical-rectangular cavities. We note that these replication rules 
are not unique and are implementation dependent. The resulting replication 
rules used in this project are given as Fortran code in Appendix C. 
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5 Scattering and Radiation 


Once the cavity aperture and volume electric fields have been determined for 
either an external excitation (scattering) or an internal excitation (radiation), 
several engineering quantities may be calculated. The aperture fields may- 
be used to determine the Radar Cross Section (RCS) for scattering or the 
radiated fields for antennas. This entails the convolution of the aperture fields 
with an appropriate Green's function. In addition, the input impedance may 
be calculated by using the interior cavity fields. In this section we will present 
the relevant formula for calculating the far-zone radiated/scattered fields and 
the input impedance from the electric fields. 


5.1 Far- field Evaluation 

Two of the most important applications of the presented formulation deal 
with the calculation of the cavity’s RCS and the radiation pattern due to 
sources placed within or on the aperture of the cavity. 

We begin with the integral representation of the scattered magnetic e 

in terms of the aperture fields. We have, 

H a (r,0,<j)) = jY 0 k o aj s G 2 {rJ,<i>-,a,<i>,z)- 

\jp{a,4> ,z) x E{a,<t>' ,z )\d<t> dz (32) 

with (r ,9,4>) indicating the observation point in spherical coordinates.^ When 
the observation point is very far from the cylinder, the dyadic Green s func- 
tion in (32) can be replaced by its far-zone representation 


(? 2 ( r, 6 , (p\ a, 4> 


i — \g h Q4^ + G dz 9z + 


(33) 


where the unprimed unit vectors are functions of the observation position an 
primed ones are functions of the integration point in (32). The components of 
this far-zone Green’s function are determined by a mode matching procedure 

giving 
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As one might expect, these series converge rather slowly for large k o as\n0. 
They must therefore be recast in another form by employing Watson's trans- 
formation and Fock theory as described before. We have. 
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where the appropriate Fock functions are given by (24). As expected, when 
d> ~ 0, the formulas attributed to Goriainov’ [15] 
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( 36 ) 


are more useful. The far- zone scattered or radiated field can be computed 
numerically by using either (34), (35), or (36) in (32). 
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For the scattering problem, the RCS is most often the quant it v of interest 
This is given by 


cr(9 , o) 


lim 

r— -oc 


' 2 \fl a (r,0.Q)\ 

\fl‘(r.O.0)\ 


(37) 


Alternatively, the antenna gain may be computed from the fa r- zone fields as 


GdB(0- <p) — 10(°6io 


4 '(i) ,|£,(, ' ,,)|, ] + I0,o * ,o & 


(38) 


where A cm is the wavelength in centimeters, 7? m is the input resistance which 
is given in the next section and E r is the radiated electric field as r — ♦ oo. 
For comparison with other techniques which define the RCS in terms of the 
electric field instead of the magnetic field in (37), these two fields may be 
related by 


E 0 = -ZJl 0 

E e = Z 0 H<, (39) 


in the far-zone. This must be considered when referring to the polarization 
of an incident or scattered field. 


5.2 Input Impedance 

In addition to the antenna pattern and gain, designers are concerned with 
the input impedance of an antenna for feedline matching purposes. The 
FEM approach allows the calculation of the input impedance of the radiating 
structure in a rather eloquent manner. The input impedance is comprised of 
two contributions [19] 

Zin = Zp + Z d (40) 

where the first term is the probe’s self-impedance which is the impedance ( 
e.g. the probe’s impedance in the absence of the patch) and the second term 
is the contribution of the patch current to the total input impedance. The 
probe self-impedance accounts for the finite radius of the probe and hence is 
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omitted when a zero-thickness probe is assumed. Ignoring the probe-feed's 
self impedance, we have [19] 


n - ~~j 2 ( E{p<0.z) • J l i ' ,t {p.o.z)pdpdodz (41) 

*o J 1 1 

where the impressed current is given by (28), \\ refers to the volume elements 
containing the probe-feed, the electric field is the interior field associated with 
the feed edge and I 0 is the constant current impressed on the probe. Utilizing 
(4) and (28) into (41) yields 
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E(>) f iPl 
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(42) 


which must be summed over the four radial edges of the element which con- 
tains the feed. 


6 Validation 


Having developed the FEM-B1 formulation for cavities recessed in an infi- 
nite, metallic cylinder and having implemented the technique in a manner 
which has low computational and storage demand, a essential task is the val- 
idation of the written computer code. Unfortunately, although their is much 
interest in the scattering and radiation characteristics of conformal patch 
arrays, a survey of the literature indicates a dearth of published data. VVe 
are currently awaiting measured data pertaining to a four-patch wraparound 
antenna mounted on a low observable test structure which is illustrated in 
figure 5. Until that data is available, we must resort to limited validation 
which either considers quasi-planar configurations or curve patches excited 
by a normally incident plane wave. We shall now look at such preliminary 
validation of the code. 

6.1 Scattering 

The first validation effort for scattering by cavity-backed patch antennas 
relies on the fact that a small patch on a very large radius cylinder is quasi- 
planar and approximates rather well an equal sized planar patch. For our 
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Unwrapped Element 

Figure 5: Wraparound test body. 


2 ’ 





t(>st we chose as a reference a planar 1.-1-1S" x 1.083" patch residing on a 
2.89 x 2.10 x 0.057 cavity filled with a dielectric having i r = 4. The 
equivalent patch on a 10A cylinder is G . 4 G° x 1.083” residing on a 12. 90° x 
2.10 x 0.057 cavity. Figure G compares this work with the planar result 
computed using a similar FEM-BI code for cavity-backed antennas recessed 
in a groundplane. Although figure G illustrates only monostatic scattering in 
the 0 = 90“ plane, additional runs for normally incident monostatic scattering 
and various bistatic situations yield similar agreement. 

Comparisons may also be made for elongated cavities and 2-D MoM re- 


sults. Fong narrow cavities have very little axial interaction for principal 
plane (0 = 90“) excitation and therefore results based on this formulation 
should compare well with corresponding 2-D data. It is well known, that the 
RCS of the 3-D scatterer of length L > A 0 is related to the corresponding 
2-D scattering of the same cross section via the relation 


03 d 



(43) 


Such a comparison is shown in figure 7 for monostatic scattering by a 45° x 
5A x 0.1 A cavity for both principal polarizations. Once again the agreement 
between the two results is excellent, thus providing a partial validation of the 
new code for highly curved geometries. We remark that similar agreement 
has been observed for bistatic scattering in the 0 = 90° plane. 

The planar approximation eliminates the effects of curvature, which is 
a primary interest in this work, and the 2-D comparisons done above are 
only valid for normal incidence. To consider oblique incidence on a highly 
curved structure, we resort to comparisons with a Body of Revolution (BOR) 
code for wraparound cavities. Since the BOR code can only model finite 
structures, we simulate an infinite cylinder by coherently subtracting the far- 
zone fields of the finite structure without a cavity from similar data which 
includes the cavity. This process is illustrated in figure 8. This procedure is 
suitable for near normal incidence and was found acceptable for near grazing 
incidence as well in the case of H-poIarization (o = 90°) case shown in figure 
9 where the data is taken for bistatic scattering in the 0 = 0° plane due to 
plane wave incidence at 0, = 90° and <f> t = 0°. 

Although the agreements with the results presented above give us con- 
fidence that our implementation yields accurate data, agreement with mea- 
sured data of the geometry shown in figure 5 would provide a more secure 


28 



Monostatic RCS (o/X ) [dBj 



Angle (0) [deg] 


Figure 6: Comparison of a planar patch (1.488 x 1.083 ) residing on a 
2.89” x 2.10" x 0.059” cavity which is filled with t T = 4 dielectric and a 
quasi-planar patch on a 10Ao cylinder. 


29 




iUMdllt {G/k ) [aK| 



Angle (0) [deg] 


Figure?: Comparison of 2-D MoM results and this work for a 45° x5A 0 x 0.1 A, 
air-filled cavity. 
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Figure 8: Far-zone subtraction scheme for the simulation of an infinite cylin 
der. 
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Bistatic RCS (c/X 2 ) [dB] 
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Observation Angle (0) [deg] 

Figure 9: Comparison of this work and processed BOR data for a 3A 0 x 0.1 A 0 
aii -filled wraparound cavity excited by a normally incident H-polarized (a = 
90°) plane wave. 




validation. The test body is a twelve inch long metallic cylinder with a radius 
of six inches. The cylinder is terminated with 53° ogival end-caps which min- 
imizes the scattering by the terminations (tips). The wraparound cavity is 
2" x 0.16" and four identical metallic patches are symmetrically place around 
the cylinder where each patch is 1.04 1 x 0.69 . 


6.2 Radiation 

As was the case for scattering, the amount of published data suitable for 
validating our formulation and associated code is quite scarce. There are 
various approximate cavity models available for coated cylinders and even 
some measured data [16, 17, IS]. However, there is no cavity-backed patch 
antenna data yet available. Nonetheless, we may compare our results with 

planar results for an initial validation. 

Figure 10 compares the co-polarized antenna pattern for a 4cm x 3cm 
patch antenna in a 6cm x 5cm x 0.0795cm cavity recessed in a metallic cylin- 
der whose radius is allowed to vary. The substrate’s permittivity is c r - 2.32 
and the feed point is at {a<p. = 2cm, r, = 1cm) relative to the center of the 
patch The small glitches in the curves associated with the small radius cylin- 
der is due to the far- zone formulas (34) and not due to the FEM formulation. 
The series (34) rely on delicate cancellations which are difficult to achieve 
numerically. Nevertheless, the agreement between the planar result and the 
case where the cylinder of radius 100 cm (« 10A at 3.17 GHz) is excellent. 
Clearly curvature has no effect on the forward direction {<j> = 0°) but the 
curvature effect is pronounced at 0 = 90°. We believe that this broadening 
of the pattern as the radius decreased may be explained by considering the 
apparent size of the patch when observed at (p = O°,0 = 90°). As the radius 
of the cylinder decreases, the patch’s planar projection becomes smaller re- 
sulting in a broader pattern. Another concern for an antenna designer is the 
effect of the cavity termination on the radiation patterns. Figure 11 shows 
the gain pattern for the same patch placed on cavities whose azimuthal ex- 
tend varies from 65° to 360° (wraparound). Obviously, near the resonant 
frequency, the effect of the cavity termination edges is only apparent near 
0 = 180° provided the cavity is small. For this test, the patch occupies 45 
of the cavity and hence the 65° cavity has only 7.5° (or 6.5 mm) to spare 
on either side of the patch. Larger cavities are associated with higher ^m 
in the rear direction (0 = ISO") possibly due to creeping wave effects. The 
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Gain [dB] 



Figure 10: Co-polarized antenna patterns for a 4cm x 3cm cm patch antenna 
on a 6cm x 5cm x 0.0795cm cavity fed at {a<t> s = 2cm, 2 = 1cm) relative to 
the center of the patch. The cavity is filled with c r = 2.32 dielectric and the 
radius of the cylinder is allowed to vary. 
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Figure 11: Effect of cavity termination on 
by the antenna used in Figure 10. 


the radiation pattern generated 




bac kvvard gain is of interest to low observable' antenna designs and we note 
that planar approximations are of no use for such antenna designs since there 
is no backward radiation. 

In addition to radiation patterns, designers are interested in input impedance 
calculations. Figure 12 compares our result for a quasi-planar patch and vali- 
dated data for the planar patch antenna considered previously. Also shown is 
the input impedance calculated as a function of frequency for a highly curved 
cylinder (a = 5 cm). The decrease in the input resistance is typical of curved 
patches (see for example [16]) and we note that the resonant frequency has 
not moved consistent with T M 0 \ mode excitation. We are currently pursuing 
the validation of our code with reference measured data for curved patches. 

7 Future Work 

In this report, we have presented a FEM-BI formulation suitable for cavity- 
backed antennas which are recessed in a circular, metallic cylinder. Along 
with the formulation, we presented an implementation strategy which mini- 
mizes both computational and storage demand. Key to this goal was the use 
of asymptotic formulas for the dyadic Green's function and the use of the 
FFT-based matrix- vector product in the BiCG solver. The resulting com- 
puter program has been validated, to the extent possible, for both scattering 
and radiation problems. 

In the near future, we want pursue further validation of the written FEM- 
BI code by comparison with measured data for curved cavity-backed anten- 
nas. This will be the first such comparison available to the electromagnetics 
community and hence will form an important resource for future work. Upon 
establishing full confidence in our implementation, we will then be able to 
study effects due to patch and cylinder of curvature on antenna parameter 
and scattering properties. We believe that the numerical laboratory de- 
veloped herein will provide a powerful tool for the analysis and design of 
conformal antenna arrays on curved surfaces. 

Having produced a code suitable for antennas mounted in metallic sur- 
faces, we want to continue this study for coated structures. The FEM-BI 
method will allow accurate computation of the input impedance, radiation 
patterns, scattering for large complex antennas mounted on inhomogeneous 
substrates. Another possible implementation could use Absorbing Bound- 
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4 cm x 3 cm Patch Antenna 



I — • — Planar 

a = 100 cm 
--*r- a = 5 cm 

Figure 12: Input impedance for the patch considered in Figure 10. 
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ary Conditions (ABC's) for mesh termination thus retaining a sparse system. 
However, the accuracy of FEM-ABC formulations for input impedance cal- 
culations have to he established and this will be an important part of this 
investigat ion. 
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Appendices 


39 




A BiCG Algorithm 

I he biconjugate gradient (BiCG) algorithm is one member of a family of iter- 
ative solvers which have proven useful in computational electromagnetics [20]. 
I he BiCG unlike the conjugate gradient (CC!) method does not guarantee 
convergence but does have the advantage of utilizing only one matrix- vector 
product in its symmetric implementation. Although the convergence char- 
acteristics of the BiCG algorithm is erratic (see for example figure A-l), it 
often conveiges in fewer iterations than the CG algorithm. I his appendix 
lists the BiCG algorithm appropriate for use with symmetric matrices [21]. 
Consider the system 

MM = M (A-i) 

wh-re (.4] is a symmetric matrix, {.r} is the unknown data vector and {6} 
is the excitation data vector. The BiCG pseudocode follows (assuming an 
initial guess {.r}, which may be {0}: 


Initialize: 

<>'), = W-MM, 
M, = W-MM, 

Iterate: 

a = < > 

< [4]{p}„.{p}; > 

t X }»+l = {*}„ + an {;>}„ 

M„+l = i r } n ~ fl n [-4] {?}„ 


__ £ )n + l ’ { 7 '}n+l 

<M„, {rj;> 


{p} — { r )n+i + C n {p} n 

(A-2) 

here the Euclidean norm is given by 


M 


< {x} ,{b} > = #[771] 6* [771] 

m= 1 

(A-3) 
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Figure A-l: Convergence behavior of a typical system solve using the BiCG 
algorithm. 
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In ( A - 2 ) the subscripts refer to the iteration and the symbol denotes 
(oniplex conjugation. Many termination criteria have been used in the past. 
One of tin’ most popular is 


< Wn + l -{rln+l > < 

<{*}•{*}> ~ f 


( A-4) 


vvheie c ^ 0.01 is a typical acceptable tolerance. As can be expected, the 
tolerance may need to be tightened or relaxed depending on the problem at 
hand and the desired accuracy of the result. Note that t should be kept small 
foi antenna impedance computations but can be relaxed tor scattering and 
radiation pattern calculations. 
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B FFT-based Matrix- Vector Product 

I lie numerical solution of integral equations (IE) (or boundary integrals) 
is often performed by converting the IE into a linear system of equations 
using the method of moments (MoM) procedure. The MoM solution often 
requires the generation and storage of 0(N' 2 ) matrix entries and, if a direct 
solver such as LU decomposition is utilized, 0( A -3 ) operations are required, 
where N is the number of unknowns. If an iterative solver such as the 
biconjugate gradient (BiCG) method is used, the solution can be found in 
0(r • / • .V 2 ) operations where r is the number of right-hand sides and i is 
the number of iterations per right-hand side. However, for eirculant or block 
circulant matrices, the solution may be reached in Q(r-i ■ .Ylog.Y) operations. 
This requires the use of FFTs for computing the matrix-vector products in 
the BiCG algorithm and consequently the resulting solution scheme is often 
leleired to as the BiCG-FFT method. In this appendix, we will present 
specific examples of this efficient technique for 2-D and 3-D geometries. We 
will first look at the relatively straightforward 2-D problem followed bv the 
necessarily more involved 3-D case. 

B.l 2-D Integral Equations 

Suppose a flat resistive strip centered at the origin of the y=0 plane is excited 
by an E-polarized [TM Z ) plane wave as shown in figure B-l. The appropriate 
Electric Field Integral Equation (EFIE) may be formed by enforcing the 
resistive transition condition on the strip giving 

E‘.(x) = R(x)J : (x) + - f J,(x )H { 0 2) (fc 0 |.r - x’|) dx (B-l) 

2 

where E z (x) is the incident electric field, R(x) is the normalized resistivity as 
a function of lateral position. J.(x) is the equivalent electric current on the 
stiip induced by the incident field and k Q — ^ is the free-space wavenumber. 
The time dependency e JW ‘ is assumed and suppressed. In (B-l), the primed 
coordinate denote the source point while the unprimed ones indicate the test 
point. Once the current has been determined by solving (B-l), the scattered 
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Figure B— 1 : Flat resistive strip geometry. 
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field in the far-zone is given by 


K(O) 



(B-2) 


where Z Q is the free-space intrinsic impedance. 

We proceed with the numerical solution of ( B— 1 ) by expanding the un- 
known current in terms of subdomain basis functions as 

JA.V) = X\7[n]ir n (.r) ( B— 3) 

n=0 


W 


/hen 


w 


u(^) = H’(.r) nA.r - — < x < (n + l)Ax - — 

= 0 else ( B— 4 ) 


and Ax = ^ 7 . The weight functions li f (x) may assume various forms, one of 
the simplest being a pulse (H 7 (x) =1). Substituting (B-3) into (B-l) and 
performing GalerkiiTs testing, we obtain 
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which is the discrete form of the integral equation. When pulse basis func- 
tions are used, (B-5) becomes 


r(n+l)Ax-f # r , /(n+l)Ax-S 

/ E\{x)dx — T J[n ]< <$[?i — n] f R(x)dx + 

JnAx— y , l JnAx—^j 
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(B-6) 
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whore the Kronecker delta function 


S[n — n } 


1 n = n 

0 n ^ n 


(B-7) 


has been introduced. Assuming that the resistivity is constant within a 
segment (e.g. R[n] = R(x) for nAx < x < (n + l)Ax) and making the 
change of variables £ = x + = x + we have 


j^ + ] 1 Et di = ^ “ "I + 

~ 71 =0 


jk /l«+l)Ai- fln' + Uii , , -I 

t/ /. 

(B-8) 


We now observe that the double integral is in convolutional form and since 
the segments are of uniform length (A.r), we may introduce the discrete 
function 


/• r(n+ 1 }Ax y(n + l)Ax . , 

</[»-»] = -f / [, II { J ] {k 0 \x - X |)d£ di (B-9) 

4 J nAx J n Ax 


and rewrite (B-8) as 

<n+l)Ax 


N - 1 


jV-1 


/ E-(0d£ = Ax X] ^[ n ]/2[«]^[n - n ] + 53 ^[ n “ n ] 

■ / " Al * n' =0 n’=0 

(B— 10) 


The first sum in (B-10) is recognized as the product of a diagonal matrix 
( Ax/?[n]<5[t7 — n'J) and the unknown data vector (J[n'],n' = 0, 1,2, ...,N — 1) 
where as the second sum is a truncated.discrete, linear convolution. 

The discrete form of the IE (B-10) may be written as a matrix equation 

\Z\{J) = {/} (B-ll) 


where the excitation vector entries are given by 

|*(n+l)Ar 
/ nAi 


/•(n + l)Ax / uj \ 

™ - L £ - : («- 2 )^ 


(B-12) 
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and l lie impedance matrix is given by 

Z[n.v] = Ax/?[n]<$[r? — n ] + (jf[n — ?; ] (B-1,'3) 

Alternatively, in preparing to take advantage of the convolution property, 
we may rewrite (13-10) as 

[Ax /?[/)]% - n ]] { J } + [ 5 [jj - »']] { J } = {/} (B-14) 

Obviously, the first matrix-vector product can be trivially computed in N 
operat ions 


//[;>] = J p [n]A.r/?[u] n = 0, 1.2 .V - 1 ( B-15) 

but the second matrix-vector product involves a fully populated matrix and 
would thus normally require 0(X 2 ) operations for its execution. However, 
since g[n - »'] is a discrete convolution operator, the product may be com- 
puted in only O(.YIogA') operations by invoking the discrete convolution 
theorem 

{.P} [<y(n - n']J = {? D {j*} • F d {g}} (B-16) 

where • denotes a Hadamard product and the discrete Fourier transform 
(DFT) pair is gi ven by 

?d {/[”»]} = fl m } e ~ J ^ rnq 

m—Q 

•^d 1 {^[7]} = j? 53 ( B— 1 7) 

7=0 

The pulse operator //[A r ] indicates that only .V of a possible 2A r - 1 values of 
the discrete convolution are retained. Also, the DFT (and its inverse) operate 
on periodic sequences of period M (all sequences with a tilde are periodic in 
this appendix). In (B-16), the DFTs must be of length M > 2N - 1 since 
a minimum of 2 N — 1 entries are needed for the complete specification of 
g[n — n ] as n and n vary from 0 to N — 1. The unknown iterate ( J p ) is 
given by 

>[n] = J[n] ii = 0, 1,2, ..., Af — 1 

0 n = N, N + 1,.V + 2 M - 1 (B-18) 
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whereas the Green’s function sequence is given 


bv 


M — 1 (B-19) 


tj[n] = </[”] n- 0.1.2 N 1 ^ 

_ q n = A, A + 1. A/ + 2, ..., 

M M , M , 0 
= <?[A/-n] n = y,y + Fy + - 

and we have assumed that </[n - n] = g[n - «]■ Note ^ at in (B " 19) ’ '[ 
M > 2N-\ it is necessary to pad g[n) with zeros in the middle of the sequence 
as shown. Also, only the first row(column) ol the Green s function matrix 
( L[„ _ n ]1 ) need be explicitly computed since this is a Toephtz matrix t us 
ensuring low 0{S) memory demand. In generating the Green's function 
sequence, tvpically the first segment (left most) is used as a source while 
the non-periodic sequence (</[»]) computed by testing at all segmen s (e.g 
t t 1 _ Q „ _ 0 1 V - 1). Thus all test segments are to the right ot 

the source' segment. If this sequence did not have the symmetry property 
la \„ _ n ’l = - 77 . 1 ), the interactions with test segments which are to 

the left of the source would need to be computed explicitly. Alternative y, 
if the sequence possessed anti-symmetry (</(» — n } — ~ s[ n " ’ P eno 

replication is still possible with the following modification to (B-19) 


g[n] = 9[”\ 

= 0 


n = 0,1,2 N- 1 


M 


n = N,N + 1,A/ “ 1 


M M . A/ _ \i — i 
= -g[M-n] n = y,y + Gy + 2 M 


(B-20) 


Combining (B-15) and (B-16), the matrix vector product in the iteration 
cycle of the BiCG-FFT method is efficiently computed as 

[Z){jn = + 

This computation requires 0(N + (M) log,( Ml) operations provided radix-2 
FFTs are used in (B-21) and should be compared to the standard matrix 

vector product 

[. Z}{J p } = £ 7 p [n']Z[n,n'] n = 0,1,2 N- 1 (B-22) 


n =0 
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Figure B-2 illustrates the comparison between (B -22) ( and discrete ' ™ nvolu ' 

, ion ( B-2 1 ) computed using radix-2 FFTs. Clearly, for . _ • - ( 1S 

more efficient than (B-22). A . 

Another 2-D geometry which yields systems which may be converted mto 

circulant matrices is the circular strip such as the one shown m figure B-3. 

The EFIE for a resistive circular strip is given by 


( / 

2 


d& d0 
( B— *23) 


where a is the radius of the arc and a is the subtended angle. Note that if 
the arc is closed, the o = 360*. A discrete form of (B-23) may be obtained 
by using pulse basis and uniform zoning (Ac = j,). We have 

r + "^£i(0-2)do = £j[n , |f«WAoi[..-n'l + 

Jn±0 ~ n ' =0 


n =0 

k a /■(n+BA* r(n' + \)M 


^ ( 2koa 
4 J J n \ 


sin 




which is similar to (B-S). The entries of the Green’s function matrix are 
given by 


M „ t L k 

9 l„_n] = — - V 


sin 




d<f> d<{> 

tR-251 


d<p 

(B 


Indeed (B-24) may be solved in exactly the same manner as (B 8) using 
FFTs. The only difference between the flat strip and the circu ar arc is a 
possible additional symmetry present in the Green’s function sequence, 
the arc is closed (a = 360°), we find that 

g[N-n] = g[n ] n = 1,2, 3 j - 1 (B- 26 ) 


which indicates that f + 1 entries of the sequence 
than N. Therefore, the solution of (B-24) for o 


need be computed rather 
= 360° requires roughly 
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Figure B-3: Circular arc geometry. 
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half the number of matrix entry evaluations as compared to an equal length 
flat strip, for n < 3G0 U . a similar symmetry exists to a lesser extent as 
long as o > ISO . In this case, although the Green s function sequence (B 
-•">) periodic, it is incorrect to assume that (B-24) now involves a circular 
convolution. The convolution is still a truncated, linear, discrete convolution 
and in practice the additional symmetry of (B-26) is not exploited unless 
matrix build time is excessive. 


B.2 3-D Integral Equations 

As was the case for 2-D geometries, there an' certain 3-D geometries which 
admit efficient solution methodologies by making use of the F FT- based im- 
plementation of matrix- vector products. Three of these geometries are the 
flat plate, an impedance insert in a ground plane and an impedance insert in 
an infinite metallic circular cylinder. 1 liese geometries are shown in figure B- 
-1. We shall now proceed with the formulation for the later two problems (e.g. 
the planar and cylindrical inserts) using a general coordinate system. In the 
following, the general coordinates (u.r) may be considered as (u = x, v = y) 
for the planar case and (u = ao. v = z) for the cylindrical case. 

An integral equation may be obtained by enforcing the standard impedance 
boundary condition (SIBC) 

h x h x E(a. i') = —})ZJi x fj(u,v) ( B— 27) 

where n is an outward directed unit normal and the total magnetic field is 
given by 

//(u.r) = H 9 °{u,v) + 

;U/C 2 (« — u,v - v‘) ■ [n' x £(i/,i/)j dudv‘( B-28) 

where A denotes the surface of the insert. In (B-2S), the geometrical optics 
fields H 3 °( u, r) is comprised of the incident field H‘(u,i >) and the reflected 

field H T (u,r). For the planar geometry, the second kind Green’s function is 
given by 


~pz- t / 

Ci2 \U — U , V — ) 


— 2 G 0 {u — u , v — v ) = 2 


7 + — 1 


e ~jk 0 \/(u-u' ) 2 

^‘ 2 -I 47 Ty/(u — u') 2 -F (u — r') 2 

(B-29) 
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Figure B-2 illustrates the comparison between (B-22) and discrete convolu- 

Thni (B-21 ) computed using radix-2 FFTs. Clearly, for A > *0, (B ) is 

mohe efficient than (B-22). ,. 

Another ‘>-D geometry which yields systems which may be converted into 

circuit matrices is the circular strip such as the one shown in figure B-3. 

The EFUSJ 01 ' a resistive circular strip is given by 


£(*) = 




sin 


0-0/ 



where a is the radiu\f the arc and a is the subtended jfrigle. Note that if 
the arc is closed, the oW 360*. A discrete form of (B-^3) may be obtained 
by using pulse basis and ihdform zoning (A<? = §)• y e have 


/ <U+ E\{0 - % )do = W[n*] fl[n]Ao% Jn ] + 
Jn±0 - n ’ -0 \ / 


k o a r( r(" + 
4 J nA<j> \ Jtx A 


Hi 2) ( 2 k 0 a 


sin 




/ N 1 
g[N-ny= $r[n] n = 1,2,3 — - 1 


(10 (10 
(B-24) 


which is similar to (B-S). The entries o^the Green’s function matrix are 
given by 

.i-.i - 

7 x (B-25) 

Indeed, (B-24) may be solved in/bxactly the same maimer as (B-8) using 
FFTs. ’ The only difference between the flat strip and the circular arc is a 
possible additional symmetry Resent in the Green’s function sequence, 
the arc is closed (a = 360°), /e find that 

\ 

(B-26) 


which indicates that f / 1 entries of the sequence need be computed rather 
than N. Therefore, yhe solution of (B-24) for o = 360° requires Whly 


50 




Figure B-4: 3-D geometries: (a) flat impedance insert in a groundplane (or 
impedance plate in free-space) (b) impedance insert in a circular cylinder 
(2a = Ad>). 
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where^ is the free-space dyadic Green's function and / = x.r + f/y + :: is the 
idem factor. For a cylinder, the appropriate Green's function is available as 
an eigenvalue series [S] for small radius cylinders or as a creeping wave series 
[!)] for large radius cylinders. Irrespective of the form of G>. upon making 
use of (B-28) into (B-2T), we have 


h x h x E(u,v) 

-Z a n x H a °(u.v) = + 


jlcji x [ G 3 (u - u'. u-v')- [n x £(«'. v')\ du dv 

(B-30) 


We note that for planar geometries, an appropriate right-handed system is 
(li.t’.n) whereas for cylindrical problems, the right-handed system is given 
by (n.u.vl Upon expanding the fields (E.H"°) and the dyadic Green's 
function (G 2 ) in terms of tangential components (B-30) yields the 

following coupled set of integral equations 


u : Z„//f(u.t’) = - 


£ u (u, v) 


+ 


v:-Z 0 H 3 u °(n,v) = - 


jk 0 J s [E.(u',i'Y C“(u-u',v-v)- 

E u (n ,l')G uv {u — u\v — r ) 
E u {u,v) 


du dv 


+ 


jkoJ[Eu{n\v)G vu {u-u\v-c)- 

E v (u\v')G vv {u - u, v - v )l du di{B-3l) 


To convert (B-31) into a discrete set of equations, we expand the unknown 
electric field components in terms of subdomain basis functions 


N v - 3 Nu-2 

E u (u, v) = £ £ E u [t,s}\V t u a {u,v) 

t = 0 3=0 

ZV„ — 2 Nu-3 

£,(«,«) = II £„MO?(u,>’) (B-32) 

(=0 3=0 
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where A u denotes the number of nodes in the (/-direction and A',, is the 
number of nodes in the e-direction. As shown in figure B-5. / denotes the 
row number of the edge discretization and s is the column number. For this 
example, the unknowns are associated with the free-edges and the associated 
basis functions may be represented as 

ll„ s ( u,v) = j — if local edge = 1 

(v t ~v b ) 

= -— — if local edge = 2 

(r t -v h ) 

0 cl. .t 

II L ! s ( u.e) = ~ — if local edge = 3 

(u r - ((,) 

( (/ — Ui ) 

= 7 i j local edge — -1 

0 else (B- 33 ) 

where the local edge numbering is illustrated in figure B-6. Note that a free 
edge is one that is not tangential to any metallic walls and that use of (B-33) 
requires a finite element type assembly. For the discretization shown in figure 
B-6. there are (,V U — 1 )(N V — 1 ) elements and a total of 2 N U N V — 3(A„ -f N v ) 
unknowns. Of these. (A r u — l)(A r v — 2) unknowns are associated with the 
(/-directed field while (A u — 2 )(N V — 1) unknowns represent the {(-directed 
field. 

Substituting the field expansions (B- 32 ) into the coupled integral equa- 
tions (B- 31 ) and employing the Galerkin’s testing procedure, we have 

= t E Eu[t\s]g uu [t-t\s-s'} + 

t ' =0 5 ; = 0 

N v — 2 AT m - 3 

E E E v [t\s'\g uv [t - i'..s - s'] t 

t'= 0 j'=0 

N , —3 /V —2 

= E E E«[t\A9vu{t-t\s-s} + 

t’ —0 3=0 


= 0,1,2 N v — 3 

= 0, 1,2, ..., N u — 2 
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(a) 



(b) 


Figure B-5: Rectangular patch discretization: (a) u-directed edges (b) v- 
directed edges. 
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'z 'z R[«’.*Vl« - >'•* - *1 1 - °- L - v - - 2 

t ' — 0 s' =0 

s = 0, 1.2 V u - 3 

( B— 34 ) 


where 




, dude 


u —u 
/ 

V = If 


Si 


-jh [ / W u {u ,v)\V u {u.v)Cr ,v (u - u',v - v')dudv'dudr 

J *• i ( 

% - e] ’ 


f elude 


u = u 
I 

y = i> 


+jA - 0 J J lV„(u\ t>)G l ’ u (t< - f - v)dudvdudv 

< Juv [t-t\*-s\ = jA-o / /, H;(u\ l /)Il' u («,t')G“ u («-u'.r-r' )du'dv'dudv 

«/ */ i / 

e 

F„[M = Z„ / W u {u,v)H 9 v °{u,v)dudv 
J St; 

F t .[/,.s] = -Z 0 [ W v [u,v)H*°{u,v)dudv (B-35) 

*/ .!><. 

and e refers to the test element while e denotes the source element. We 
note that each of the double sums in (B-34) is a truncated, discrete, linear 
convolution and hence amenable to the BiCG-FFT method. 

To proceed, we define the 2-D DFT pair (analogous to (B-17)) 

A/l-1 Mj-l _ 2n / , , , 

{/[<,*]} = E E /Me"'® 1 ” 5 ' 

t=0 3=0 

*55 {*M - JTXT *E *E (B-36) 

*• J M 1-1/2 (= o 3=0 

Using (B-36), the convolutions in (B-34) can then be rewritten as 

A £ 1A E E[t\s']g\t-t\s-s] = ^2D {^2D{E[t,s])9^2D{g[U»]}} 

t'= 0 s =0 
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(B-37) 


TIk> order of the relevant DFTs must he > 2(number of rows) - 1 and 
.\l 2 > -(tuiiiiher of columns) - 1 where the number of rows and columns 
of the discretization may vary with each convolution (see figure B5). For 
example, the first convolution in ( B-34 ) is associated with (‘/-testing and 
//-source edges and hence the number of rows and columns is ( ,V t , — 2) and 
( A u — 1 ). respectively. The field sequences are loaded into a M , x \/ 2 array in 
row/column order of the field discretization and the remaining entries form 
a zero pad. The Green’s function sequence must be loaded into a similar 
array (in the same manner) and periodic replication must be performed to 
provide the necessary "negative lags”. If the sequence has the property. 
g[l — t , .s — .* ] = g\t — t. s — .$], then this replication process takes the form 

//[/.. s] = //[/,. s] 0<f<^-10<s<^-l 

= g[M x + 2 — t. $} ^ < / < Mi - 1 0 < .s < - 1 

= g[t. \1 2 + 2 - s] 0 < t < ^ - 1 ^ 5 < A/ 2 - 1 

= <j[Mi +2 - t. M 2 + 2 - s] < t < Mi - l < s < M 2 - 1 

(B-38) 

For the presented example, and possess this property while 

the cross-terms, f/ U t/[C- s ] and g vu [t<$}, do not. If anti-symmetry is present 
then a more complex replication scheme similar to (B-20) may be used. 
Otherwise, all possible lags must be computed requiring longer matrix build 
time since more than the first u-directed and r-directed edges need be used 
as sources. 

Once the periodic arrays are loaded, the required matrix-vector product 
for the tiii-interactions may be performed in C((M x \ogM x )( A/ 2 logA/ 2 )) op- 
erations rather than the 0( ( ( A r u — 2 )(N V — 3)) 2 ) operations required for a 
standard matrix- vector product. The comparison is shown in figure B-7 with 
M x = 2(N V — 3), A/ 2 = 2( — 2) and N v = N u = N. Clearly, when the 
number of nodes per side exceeds 10-15, the FFT-based matrix-vector prod- 
uct is more efficient than a conventional matrix-vector product. In practice, 
the FFT-based product is more efficient than a standard product in terms 
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Complexity 



Figure B-7: Comparison of operation count lor a full matrix-vector product 
verses a FFT-based matrix-vector product. N is the number of nodes in each 
direction of the grid, Mi = 2(N-3) and M 2 = 2(N-2). 
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ol uallclock tunc for A < 10 since in order to exploit the memory savings 
afforded by uniform zoning of a convolutional kernel without using FFTs. ad- 
ditional overhead is incurred to match the appropriate matrix entry with the 
correct vector entries. Similar results are obtained for the other convolutions 
in ( B-34). 

C Replication Rule (Fortran) 

I his appendix contains the Fortran source code which performs the periodic 
replication of the four components of the boundary integral sub-matrix. It 
is included to illustrate the special care required for the cross-terms (0 — z 
and c - o) as opposed to the application of (B-3S) for the like-terms { 0-0 
and r — r). 
c 

c Augment the arrays 


c Provide negative lags for like-terms 
c 

Do row = l.rowFFT 

Do column = l.colFFT 

If( ( row. LE.rowFFT/2). AND. (column. LE.colFFT/2)) Then 

c First quadrant, do nothing ... 

Elself( (row.CT.rowFFT /2). AND. (column. LE.colFFT /2)) Then 
c Second quadrant 

g l T U 2 D ( row , col u m n ) = gUU2D(rowFFT+2 -row, col umn) 
gVV2D( row, column) = gVV2D(rowFFT+2-row, column) 
ElseIf((row. LE.rowFFT/2). AND. (column. GT.colFFT/2)) Then 

c Third quadrant 

gUU2D( row, column) = gUU2D(row.colFFT+2-column) 
gVV2D( row, column) = gVV2D(row.colFFT+2-col umn) 

Else 

c Fourth quadrant 

gl T U2D( row, column ) = 

gUU2D( row F FT +2-ro\v,coIFFT+2-column) 
gV V2D(row,column) = 
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e v V2 D ( row F F T+2- row .col F FT + 2-col u m n ) 

Encllf 

EndDo 

EndDo 

c 

c Special treatment for cross-terms 
c 

It'( .NOT. wrapAround) Then 
Do row = l.rowFFT/2 

Do column = l.colFFT/2 


c Replicate lor 1 V 

c 

c Fourth quadrant from first quadrant 

gf V2d ( row FFT+ 1 -row. col FFT+1 -column) = 

gUV2d(row, column) 
c Second quadrant from first quadrant 

gl ■ V2d( rowFFT + 1 - row. column ) = -guv2d(row, column) 

c 

c Replicate lor \ U 
c 

c Fourth c|uadrant from first quadrant 

g V U 2d ( row F FT + 1 - row.col FFT + 1 -column ) = 
g V U 2d ( row ,col u mn ) 
c Third quadrant from first quadrant 

gVU2d(row,colFFT+ 1-column) = -gVU2d(rovv, column) 

End Do 
EndDo 
c 

Do row = 1 .row FFT /2 

Do column = l,colFFT/2 
c 

c Now mirror into first and third quadrants foi UV 

c 

c First quadrant from second quadrant 

gUV2d(row, column) = (0.0D0.0.0D0) 

g UV2d( row, column) = - g UV2d(rowFFT-l-row, column) 
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c I bird quadrant from fourth quadrant 

gl A r 2d(ro\v,colFFT+l-colunin) = 

-gl V2d(ro\vFFT- l-ro\v,coIFFT + 1 -column) 

c 

c Now mirror into first and second quadrants for VU 
(' 

c Fi rst quadrant from third quadrant 

gVU2d(row, column) = (0.0D0,0.0D0) 
gVU2d(row, column) = -gVU2d( row. colFFT- 1-column ) 
c Second quadrant from fourth quadrant 

g\*l 7 2d(rowFFT + 1-row. column ) = 

-gVU2d(rowFFT + 1-row.colFFT- 1 -column) 

Enel Do 
EndDo 
Else 

c WRAPAROUND CAVITY 
Do row = 1. row F FT/ 2 

Do column = l,colFFT/2 

c 

c Replicate for UV 
c 

c Fourth quadrant from first quadrant 

g U V 2d ( row F FT+1 - row ,col F FT +1 - co 1 u m n ) = 
g U V2d ( row ,col um n ) 
c Second quadrant from first quadrant 

g U V2d ( row FFT+ 1 - row .col u m n ) = 
gU V2d( row ,guvColMax+l -column) 
c 

c Replicate for VU 
c 

c Fourth quadrant from first quadrant 

gVU2d(rowFFT+ 1 -row, col FFT+1 -column) = 
gVU2d( row, column) 
c Third quadrant from first quadrant 

gVU2d(row,colFFT-fl-column) = -gVU2d(row,column) 
EndDo 
EndDo 
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Do row = l,rowFFT/2 

Do column = l,colFFT/2 

c 

c Now minor into first and third quadrants for UV 
c 

c First quadrant from second quadrant 

gfi'V2d(row, column) = (O.ODO,O.ODO) 
g U V 2d ( row .column) = -gU V2d ( row F FT- 1 - row , col u mn) 
c Third quadrant from fourth quadrant 

gU V2d( row, col FFT+1 -column) = 

-gl \ 2d ( row F FT- 1-row.colFFT + 1 -column) 
c 
c 

c Now mirror into first and second quadrants for VU 
c 

c First quadrant from third quadrant 

g\U2d( row, column) = (0.0D0.0.0D0) 
g\’U2d( row. column) = 

gVU2d(row,colFFT-gvuColMax + column) 
c Second quadrant from fourth quadrant 

gVU2d(rowFFT+ 1-row, column) = 

gVU2d(rowFFT+Trow,colFFT-gvuColMax -(- column) 

End Do 
EndDo 
Endlf 
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